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Abstract 

An  algorithm  is  developed  based  on  the  multi-fluid  plasma  model  derived 
from  moments  of  the  Boltzmann  equation.  The  MHD  (magnetohydrodynamic) 
model  involves  assumptions  that  limit  its  applicability.  The  multi-fluid  plasma 
model  only  assumes  local  thermodynamic  equilibrium  and,  therefore,  accurately 
models  the  appropriate  physical  processes.  The  multi-fluid  plasma  model  typi¬ 
cally  has  two  fluids  representing  electron  and  ion  species.  Large  mass  differences 
between  electrons  and  ions  introduce  disparate  temporal  and  spatial  scales  and 
require  a  numerical  algorithm  with  sufficient  accuracy  to  capture  the  multiple 
scales.  The  multi-fluid  capability  is  not  limited  to  two  species.  Plasma  with 
multiple  components  can  be  modeled,  e.g.  impurity  ions,  neutral  gas.  The 
multi-fluid  equations  are  derived  in  divergence  form  for  the  naturally  occurring 
conserved  variables.  The  source  terms  of  the  multi-fluid  plasma  model  couple 
the  fluids  to  themselves  (interspecies  interactions)  and  to  the  electromagnetic 
fields.  The  solution  and  evolution  must  be  tightly  coupled  to  prevent  unstable 
numerical  oscillations.  The  electromagnetic  field  equations  are  solved  by  apply¬ 
ing  correction  potentials  or  by  using  a  mixed  potential  formulation  to  eliminate 
any  components  of  the  fields  that  violate  the  divergence  constraints.  A  dis¬ 
continuous  Galerkin  method  is  developed  to  solve  the  governing  equations  on 
a  computational  grid  and  to  simulate  plasma  phenomena.  Interspecies  inter¬ 
actions  also  occur  through  collisional  source  terms  that  account  for  the  direct 
transfer  of  momentum  and  energy.  In  addition  to  the  plasma  and  electrody¬ 
namic  physics,  the  multi-fluid  plasma  model  captures  atomic  physics  in  the 
form  reaction  rate  equations  for  ionization  and  recombination,  which  introduce 
new  temporal  scales  to  the  plasma  dynamics  model.  The  numerical  algorithm 
must  treat  the  inherent  stiffness  introduced  by  the  multiple  physical  effects  of 
the  model.  Fully  implicit  and  semi-implicit  treatments  have  been  implemented 
with  the  semi-implicit  treatment  showing  promise  for  a  more  optimum  method. 
Nonreflecting  boundary  conditions  using  a  lacunae-based  method  have  been  im¬ 
plemented  and  provide  higher  solution  fidelity  for  open  boundary  problems. 
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1  Executive  Summary 

This  project  represents  a  two-year  effort  to  develop  a  new  algorithm  for  plasma 
simulations  based  on  the  multi-fluid  plasma  model.  The  algorithm  is  capable  of 
three-dimensional,  time-dependent  plasma  simulations  that  capture  the  multi¬ 
physics  and  multi-scale  effects  of  the  complete  multi-fluid  plasma  model.  The 
multi-fluid  plasma  equations  is  formulated  in  divergence  form  which  includes 
source  terms.  The  source  terms  couple  the  fluids  to  the  electromagnetic  fields, 
account  for  collisions  between  the  fluids,  and  provide  sinks  and  sources  from  the 
atomic  reactions.  A  high-order  algorithm  is  developed  that  solves  the  complete 
multi-fluid  plasma  model  such  that  equilibrium,  perturbation,  and  transient 
phenomena  can  be  accurately  simulated.  The  algorithm  uses  the  discontinuous 
Galerkin  method  to  achieve  high-order  accuracy.  Novel  numerical  approaches 
for  time  integration  are  investigated  to  perform  simulations  that  are  not  limited 
by  the  shortest  timescale  of  the  system.  These  methods  include  fully  implicit 
backward  difference  and  a  semi- implicit  treatment.  The  algorithm  takes  ad¬ 
vantage  of  the  parallel  computing  architectures  that  are  available  locally  and 
at  the  Major  Shared  Resource  Centers  (parallel  workstation  cluster,  Dell  Pow- 
erEdge,  Cray  XD1,  and  others).  The  new  algorithm  has  been  benchmarked 
against  known  analytical  results,  previously  published  results,  and  experimen¬ 
tal  data  from  the  Air  Force  Research  Laboratories  [the  field  reversed  configu¬ 
ration  (FRC)  implosions  for  magnetized  target  fusion  (MTF)  at  Kirtland  AFB 
and  FRC  experiments  at  the  University  of  Washington]. 

Many  plasma  simulation  codes  are  currently  based  on  the  MHD  (magne¬ 
tohydrodynamic)  model.  The  derivation  of  the  MHD  model  involves  several 
assumptions  that  severely  limit  its  applicability.  MHD  simulations  of  several 
technologies  important  to  the  USAF  have  failed  to  predict  the  experimentally 
observed  plasma  behavior,  e.g.  stability  of  field  reversed  configurations.  Recent 
success  with  the  two-fluid  plasma  model  has  provided  a  development  path  for 
plasma  simulations  that  are  more  physically  accurate  and  capable  than  MHD 
models.  [1-6]  Building  on  the  successful  development  of  the  two-fluid  model, 
a  plasma  simulation  algorithm  is  developed  for  the  multi-fluid  plasma  model. 
The  multi-fluid  plasma  model  only  assumes  local  thermodynamic  equilibrium 
within  each  fluid  and  is,  therefore,  more  physically  accurate  and  capable  than 
MHD  models.  The  separate  electron  and  ion  response  is  accurately  represented, 
as  well  as,  multi-component  plasmas  that  can  include  impurities,  neutrals. 

The  multi-fluid  model  is  formulated  in  conservation  form  and  lends  itself 
naturally  to  accurate  fluid  models.  For  example,  an  approximate  Riemann 
solver  is  developed  for  the  two-fluid  plasma  model  to  compute  the  fluxes  in 
a  stable  and  accurate  manner.  Several  methods  are  investigated  to  solve  the 
electromagnetic  field  model,  which  includes  the  source  terms  and  divergence 
constraints.  These  methods  include  a  purely  hyperbolic  formulation  of  the 
Maxwell’s  equation  and  a  mixed  potential  formulation.  The  plasma  fluids  (elec¬ 
trons  and  ions)  and  the  electromagnetic  fields  communicate  through  the  source 
terms.  The  fluid  momentum  and  energy  equations  have  source  terms  that 
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depend  on  E  and  B.  The  electromagnetic  equations  have  source  terms  that 
depend  on  v,  and  ve  (Ampere’s  law)  and  nt  and  ne  (Gauss’s  equation).  All  flu¬ 
ids  (plasma  and  neutral)  communicate  through  scattering  collisions  and  atomic 
reactions  that  appear  as  source  terms.  Accurately  coupling  the  source  terms  is 
important  both  for  numerical  stability  and  for  modeling  plasmas  where  large 
equilibrium  forces  exist. 

The  multi-fluid  algorithm  uses  a  discontinuous  Galerkin,  finite  element 
method  for  the  spatial  representation  and  a  TVD  Runge-Kutta  method  for 
the  time  advance.  [2]  Solutions  are  found  with  up  to  16th  order  spatial  accu¬ 
racy  and  3rd  order  temporal  accuracy.  [7]  The  multi-fluid  plasma  algorithm 
is  used  to  model  multiscale  physics  of  current-carrying  plasmas,  such  as  the 
Z-pinch  [4]  and  the  field  reversed  configuration  (FRC)  [5].  These  plasma  con¬ 
figurations  balance  large  equilibrium  forces  between  the  plasma  pressure  and 
the  electromagnetic  pressure.  The  high-order  algorithm  is  seen  to  significantly 
improve  the  ability  to  maintain  equilibrium  with  no  artificial  decay.  [7] 

Domain  truncation  is  accomplished  by  implementing  open  boundary  con¬ 
ditions  that  are  based  on  lacunae  methods.  The  open  domain  boundary  is 
appended  with  an  exterior  domain  that  matches  the  interior  solution  and  then 
damps  the  solution  before  artificial  reflections  develop.  Multi-fluid  and  arbitrar¬ 
ily  complex  geometry  capability  has  been  demonstrated.  [8]  Plasma  formation 
process  is  modeled  in  a  two-dimensional  geometry  which  produces  a  plasma 
from  the  collisions  of  electrons  and  neutrals.  During  the  initial  formation  a 
Langmuir  wave  is  observed  to  propagate.  FRC  simulations  are  performed  in 
a  five-block  cylindrical  grid.  The  simulation  results  show  the  development  of 
short  wavelength  drift  turbulence. 

The  divergence  constraints  of  Maxwell’s  equations  can  be  difficult  to  satisfy 
with  the  presence  of  current  and  charge  sources  on  an  arbitrary  computational 
grid.  The  divergence  constraints  are  satisfied  by  reformulating  Maxwell’s  equa¬ 
tions  to  include  correction  potentials.  The  approach  involves  coupling  the  di¬ 
vergence  constraint  equations  with  the  time-dependent  field  equations  to  form 
a  purely  hyperbolic  equation  set.  [9]  An  alternative  formulation  of  Maxwell’s 
equations  using  mixed  potential  is  also  implemented.  The  mixed  potential  for¬ 
mulation  automatically  satisfies  the  divergence  constraints;  however,  a  gauge 
condition  must  then  be  enforced. 

This  project  was  performed  by  Prof.  Uri  Shumlak  and  graduate  students 
Robert  Lilly,  Noah  Reddell,  Eder  Sousa,  Bhuvana  Srinivasan,  and  Andree  Su- 
santo.  This  project  resulted  in  doctoral  dissertations  and  master  thesis: 

•  Bhuvana  Srinivasan,  “Numerical  Methods  for  3-dimensional  Magnetic 
Confinement  Configurations  using  Two-Fluid  Plasma  Equations”,  Ph.D. 

2010. 

•  Andree  Susanto,  “Development  of  Electromagnetic  Solvers  for  Use  with 
The  Two-Fluid  Plasma  Algorithm”,  M.S.  2009. 

These  dissertations  and  theses  can  be  obtained  from  the  University  of  Washing¬ 
ton  library  system  or  from  the  project  website,  http:/ /www. aa.washington.edu/research/cfdlab/. 
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Archival  journal  and  conference  papers  were  published  reporting  on  the  work 
from  this  project: 

•  U.  Shumlak,  R.  Lilly,  N.  Reddell,  E.  Sousa,  and  B.  Srinivasan,  “Advanced 
physics  calculations  using  a  multi-fluid  plasma  model” ,  Computer  Physics 
Communications  182  1767-1770  (2011). 

•  B.  Srinivasan,  A.  Hakim,  and  U.  Shumlak,  “Numerical  methods  for  two- 
fluid  dispersive  fast  MHD  phenomena”,  Communications  in  Computa¬ 
tional  Physics  10  (3),  183-215  (2011). 

•  J.  Loverich,  A.  Hakim,  and  U.  Shumlak,  “A  discontinuous  Galerkin  method 
for  ideal  two-fluid  plasma  equations” ,  Communications  in  Computational 
Physics  9,  240-268  (2011). 

•  A.  Hakim  and  U.  Shumlak,  “Two-fluid  physics  and  field-reversed  config¬ 
urations”,  Physics  of  Plasmas  14,  055911  (2007). 


2  Project  Results 

Plasmas  are  essential  to  many  technologies  that  are  important  to  the  Air  Force, 
some  of  which  have  dual-use  potential.  These  applications  include  portable 
pulsed  power  systems,  high  power  microwave  devices,  drag  reduction  for  hyper¬ 
sonic  vehicles,  advanced  plasma  thrusters  for  space  propulsion,  nuclear  weapons 
effects  simulations,  radiation  production  for  counter  proliferation,  and  fusion 
for  power  generation.  In  general,  plasmas  fall  into  a  density  regime  where  they 
exhibit  both  collective  (fluid)  behavior  and  individual  (particle)  behavior.  The 
intermediate  regime  complicates  the  computational  modeling  of  plasmas. 


2.1  Plasma  Models  —  Kinetic,  PIC,  MHD,  Multi- 
Fluid 


Plasmas  may  be  most  accurately  modeled  using  kinetic  theory.  The  plasma  is 
described  by  distribution  functions  in  physical  space,  velocity  space,  and  time, 
/(x,  v,  t).  The  evolution  of  the  plasma  is  then  modeled  by  the  Boltzmann 
equation. 


dfa 

dt 


+  v0 


%  +  (E 

ax  mn 


v«  x  B  ) 


dfa 

<9v 


dfa 

dt 


(1) 


collisions 


for  each  plasma  species  a  =  ions,  electrons.  The  Boltzmann  equation  cou¬ 
pled  with  Maxwell’s  equations  for  electromagnetic  fields  completely  describe 
the  plasma  dynamics.  [10-12]  However,  the  Boltzmann  equation  is  seven  di¬ 
mensional.  As  a  consequence  of  the  large  dimensionality  plasma  simulations 
using  the  Boltzmann  equation  are  only  used  in  very  limited  applications  with 
narrow  distributions,  small  spatial  extent,  and  short  time  durations.  [13, 14] 
The  seven  dimensional  space  is  further  exacerbated  by  the  high  velocity  space 
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that  is  unused  except  for  tail  of  the  distribution  or  energetic  beams.  Boundary 
conditions  are  difficult  to  implement  in  kinetic  simulations. 

Particle  in  cell  (PIC)  plasma  model  apply  the  Boltzmann  equation  to  repre¬ 
sentative  superparticles  which  are  far  fewer  (107)  than  the  number  of  particles 
in  the  actual  plasma  (1020).  [15]  PIC  simulations  have  similar  limitations  as 
simulations  using  kinetic  theory  due  to  statistical  errors  caused  by  the  fewer 
superparticles.  Boundary  conditions  are  also  difficult  to  implement  in  PIC 
simulations. 

The  other  end  of  the  spectrum  in  plasma  model  involves  taking  moments  of 
the  Boltzmann  equation  and  averaging  over  velocity  space  for  each  species  which 
implicitly  assumes  local  thermodynamic  equilibrium.  The  resulting  equations 
comprise  the  multi-fluid  plasma  model,  where  a  fluid  is  defined  for  each  species, 
e.g.  electrons,  ions,  neutrals.  When  only  two  fluids,  electrons  and  ions,  are 
included,  the  resulting  system  is  the  two-fluid  plasma  model.  [1]  The  two-fluid 
equations  can  be  combined  to  form  the  MHD  model.  [16]  However,  in  the 
process  several  approximations  are  made  which  limit  the  applicability  of  the 
MHD  model  to  low  frequency  and  ignores  the  electron  mass  and  finite  Larmor 
radius  effects. 

The  MHD  model  treats  the  plasma  like  a  conducting  fluid  and  assigning 
macroscopic  parameters  to  describe  its  particle-like  interactions.  Plasma  simu¬ 
lation  algorithms  based  on  the  MHD  model  have  been  very  successful  in  model¬ 
ing  plasma  dynamics  and  other  phenomena.  Codes  such  as  MACH2  are  based 
on  arbitrary  Lagrangian/Eulerian  formulations.  [17]  ALE  codes  are  well  suited 
for  simulating  plasma  phenomena  involving  moving  interfaces.  [18]  However, 
ALE  codes  cannot  be  formulated  as  conservation  laws  and  lack  many  of  the 
inherent  conservative  properties.  The  MHD  model  has  been  successfully  im¬ 
plemented  in  conservative  form  to  simulate  realistic  three-dimensional  geome¬ 
tries.  [19,20] 

A  severe  limitation  of  the  MHD  model  is  the  treatment  the  Hall  effect 
and  diamagnetic  terms.  These  terms  represent  the  separate  motions  of  the 
ions  and  electrons.  The  Hall  effect  and  diamagnetic  terms  also  account  for 
ion  current  and  the  finite  ion  Larmor  radius.  These  effects  are  important  in 
many  applications  such  as  electric  space  propulsion  thrusters:  Hall  thrusters, 
magnetoplasmadynamic  (MPD)  thrusters,  Lorentz  force  thrusters.  The  Hall 
term  is  also  believed  to  be  important  to  electrode  effects  such  as  anode  and 
cathode  fall  which  greatly  affect  many  directly  coupled  plasma  applications. 
Furthermore,  the  Hall  and  diamagnetic  effects  may  be  important  for  hypersonic 
flow  applications.  [21] 

The  Hall  terms  can  be  difficult  to  stabilize  because  they  lead  to  the  whistler 
wave  branch  of  the  dispersion  relation.  The  phase  and  group  velocities  of  the 
whistler  wave  increase  with  frequency.  The  velocities  become  large  even  for 
modest  values  of  the  Hall  parameter.  See  Fig.  1  for  the  dispersion  diagram. 

A  semi-implicit  technique  has  been  applied  to  treat  the  Hall  term  in  a  Hall- 
MHD  model.  [22,  23]  After  the  hyperbolic  terms  of  the  MHD  equations  are 
advanced,  the  Hall  terms  are  treated  independently.  The  conserved  variables 
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Figure  1:  Dispersion  relations  for  the  two- fluid  plasma  model  and  for  the  Hall-MHD 
plasma  model  that  results  when  asymptotic  approximations  are  applied  to  the  two- 
fluid  plasma  model.  For  small  wave  numbers  and  low  frequencies  (right  plot),  the 
upper  branch  of  the  Hall-MHD  wave  follows  the  R  wave  of  the  two-fluid  model. 
However,  the  waves  diverge  and  the  Hall-MHD  wave  fails  to  follow  the  resonance 
at  the  cyclotron  frequency.  The  wave  speed  grows  without  bound.  Artificial  hyper¬ 
resistivity  is  required  to  damp  this  branch  of  the  Hall-MHD  wave. 
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are  then  corrected.  The  procedure  can  be  computationally  intensive.  The  oper¬ 
ator  stencil  uses  5  points  in  the  sweep  direction  and  3  points  in  each  orthogonal 
direction.  The  complete  operator  stencil  is  45  points.  The  semi-implicit  method 
works  adequately  for  small  Hall  parameters,  but  becomes  unstable  or  slow  to 
converge  for  the  large  Hall  parameters  often  seen  in  applications. 

As  mentioned  above,  the  multi-fluid  plasma  model  is  more  complete  than 
either  the  MHD  or  Hall-MHD  model.  The  multi-fluid  plasma  model  resolves 
plasma  oscillations  and  speed  of  light  propagation.  However,  many  applications 
are  adequately  modeled  by  lower  frequency  dynamics.  Asymptotic  approxima¬ 
tions  (me  — >  0,  c  — >  oo)  have  been  applied  to  the  two-fluid  plasma  model  to 
eliminate  the  high  frequency  waves  that  limit  the  maximum  numerical  time 
step.  Neglecting  electron  inertia  removes  the  limitation  due  to  the  electron 
plasma  and  cyclotron  frequencies.  Infinite  light  speed  removes  the  limitation 
due  to  light  transit  times.  The  asymptotic  approximations  reduce  the  two-fluid 
plasma  model  to  the  Hall-MHD  model.  However,  applying  these  approxima¬ 
tions  fundamentally  changes  the  dispersion  relation,  as  evident  in  Fig.  1,  and 
introduces  unphysical  wave  behavior.  Specifically,  the  phase  and  group  veloc¬ 
ities  of  a  Hall-MHD  wave  increase  without  bound  with  wave  number.  The 
large  wave  speeds  increases  the  stiffness  of  the  equation  system  making  accu¬ 
rate  numerical  solutions  difficult.  Furthermore,  the  maximum  wave  number  is 
usually  set  by  either  the  computational  mesh  spacing  ( kmax  oc  Ax)  or  by  an 
artificial  resistivity.  Rigorous  convergence  studies  are  difficult  with  the  simpler 
plasma  models  since  decreasing  Ax  leads  to  larger  kmax  and  shorter  wavelength 
phenomena. 

2.2  Multi-Fluid  Plasma  Algorithm 

The  complexity  of  the  multi-fluid  model  is  greater  the  MHD  model  but  signif¬ 
icantly  less  than  the  kinetic  model.  In  this  project  a  new  algorithm  is  devel¬ 
oped  that  solves  the  multi-fluid  plasma  model  using  an  approximate  Riemann 
solver.  [1]  The  method  tracks  the  wave  propagation  across  the  domain  based 
on  conservation  laws. 

The  governing  equations  for  the  multi-fluid  plasma  model  (including  the 
electromagnetic  equations)  are  expressed  in  divergence  form.  Similar  to  the 
two-fluid  plasma  model,  the  governing  equations  for  the  multi-fluid  plasma 
models  are  derived  by  taking  moments  of  the  Boltzmann  equation,  Eq.  (1),  for 
each  species.  Local  thermodynamic  equilibrium  within  each  fluid  is  assumed. 
The  fluid  variables  are  derived  by  taking  moments  of  the  distribution  function. 

The  evolution  of  the  particle  densities  is  expressed  by  continuity  equations. 
The  equations  are  the  zeroth  moment  of  the  Boltzmann  equation. 

dna  ( FA  _  dna 

dt  \Qa)  dt 

where  na  is  the  number  density  for  species  a  and  the  particle  fluxes  are  defined 
by  the  partial  current  densities  jQ  =  <7anava  in  terms  of  the  charge  qa  and  fluid 


(2) 
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velocity  va  for  each  species.  The  net  production  rate  of  species  a  due  to  atomic 
reactions  is  denoted  on  the  right-hand  side  of  the  equation  with  a  subscript  T. 
Contributions  due  to  atomic  reactions  are  described  later. 

The  first  moment  of  the  Boltzmann  equation  yields  momentum  equations 
for  each  species.  The  momentum  equations  are  written  in  divergence  form  in 
terms  of  the  partial  current  densities. 
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where  E  and  B  are  the  electric  and  magnetic  fields,  PQ  =  pa I  +  na  is  the  total 
pressure  tensor  of  species  a  (sum  of  the  scalar  pressure  and  the  stress  tensor), 
and  Ka0  is  the  momentum  transfer  vector  from  species  a  to  species  (3.  The 
net  current  density  (momentum)  generation  rate  of  species  a  due  to  atomic 
reactions  is  denoted  on  the  right-hand  side  of  the  equation  with  a  subscript  T. 

The  second  moment  of  the  Boltzmann  equation  yields  energy  equations  for 
each  species  which  are  expressed  in  divergence  form  for  the  total  energy. 
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where  the  total  energy  is  defined  by 


£a  = - TPa  +  -manaVa 

7  —  i  z 


(4) 

(5) 


where  7  is  the  ratio  of  specific  heats.  An  adiabatic  equation  of  state  is  assumed. 
The  energy  addition  rate  of  species  a  due  to  atomic  reactions  is  denoted  on  the 
right-hand  side  of  the  equation  with  a  subscript  T  and  is  described  in  a  later 
section. 

Neutral  fluids  can  be  incorporated  into  the  multi-fluid  plasma  model  by 
eliminating  the  species  charge  which  returns  the  conventional  expressions  for 
the  governing  equations. 


2.2.1  Electromagnetic  Field  Evolution 


The  electromagnetic  fields  influence  the  motion  of  the  plasma  fluid  through  the 
Lorentz  force  which  is  contained  in  Eq.  (3).  The  motion  of  the  plasma  influ¬ 
ences  the  evolution  of  the  electromagnetic  fields  through  the  redistribution  of 
charge  density  and  current  density.  Maxwell’s  equations  govern  the  evolution 
of  the  electromagnetic  fields.  The  net  charge  density  Yla  Qa^a  &nd  total  cur¬ 
rent  density  ja  are  calculated  directly  from  the  multi-fluid  equations  which 
couple  the  electromagnetic  fields. 
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e0V  •  E  =  ^  qanQ 

Ot 

(8) 

V-B  =  0 

(9) 

An  important  area  of  progress  is  investigating  accurate  electromagnetic  field 
solvers.  The  electromagnetic  field  model  includes  divergence  constraint  rela¬ 
tions,  which  if  not  accurately  satisfied,  can  lead  to  nonphysical  solutions.  Spe¬ 
cial  treatment  is  required  because  the  divergence  relations  over-constrain  the 
solution.  Satisfying  the  divergence  constraint  relations  requires  adding  correc¬ 
tion  potentials  to  form  purely  hyperbolic  equations  [9],  which  requires  solving 
additional  hyperbolic  equations  to  sweep  the  divergence  error  out  of  the  domain 
or  formulating  Maxwell’s  equations  with  mixed  potentials  (scalar  and  vector 
potentials)  that  automatically  satisfy  the  divergence  constraint  relations. 

The  purely  hyperbolic  version  of  Maxwell’s  equations  are  expressed  as 


C/J-> 

— — b  V  X  E  +  7V i/j  =  0 
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(10) 

c2V  X  B  +  xc2V<)>  = - qanava 
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(12) 

^+7c2V-B  =  0 
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(13) 

where  error  propagation  speeds  7  and  x  are  introduced  to  convect  the  diver¬ 
gence  errors  out  of  the  domain. 

The  electromagnetic  fields  can  be  formulated  using  potentials  that  mathe¬ 
matically  satisfies  the  divergence  constraint  conditions,  Eqs.  (8)  and  (9).  The 
mixed  potential  formulation  is  derived  by  expressing  the  electric  and  magnetic 
fields  using  scalar  and  vector  potentials.  Specifically,  the  scalar  potential  <i>  and 
vector  potential  A  are  defined  by 


^  dA 

E  =  -v$-  «-• 

(14) 

B  =  V  x  A. 

(15) 

The  potential  formulation  requires  setting  an  arbitrary  gauge  condition  that 
affects  the  form  of  the  resulting  governing  equations.  If  a  Coulomb  gauge 
condition  (V  •  A  =  0)  is  assumed,  the  evolution  equations  become 
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Figure  2:  Evolution  of  the  magnetic  field  divergence  for  the  purely  hyperbolic  and 
mixed  potential  formulation  of  Maxwell’s  equations. 


If  the  Lorenz  gauge  V  ■  A  =  —  is  assumed,  the  field  equations  become 
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Both  the  purely  hyperbolic  and  the  mixed  potential  formulations  have  been 
implemented  to  evaluate  the  ability  of  each  method  to  preserve  the  diver¬ 
gence  constraint  relations.  Figure  2  shows  the  results  for  the  GEM  challenge 
(Geospace  Environmental  Modeling  Magnetic  Reconnection  Challenge)  prob¬ 
lem  of  collisionless  magnetic  reconnection.  [24]  The  mixed  potential  formula¬ 
tions  show  a  divergence  error  on  the  order  of  machine  precision  and  much  lower 
than  the  purely  hyperbolic  formulation.  However,  the  computational  effort  for 
the  mixed  potential  formulation  is  significantly  greater. 


2.2.2  Characteristics  of  the  Multi-Fluid  Plasma  Model 

As  presented  above,  the  governing  equations  for  the  multi-fluid  plasma  model 
including  the  electromagnetic  equations  are  expressed  in  divergence  form. 

^  +  V  •  F  =  5  (20) 

where  Q  is  the  vector  of  conserved  fluid  variables,  F  is  the  tensor  of  hyperbolic 
fluxes  and  S  is  the  vector  containing  the  source  terms.  The  length  of  vectors 
depends  on  the  number  of  fluids  included  in  the  model.  For  Q,  there  are 


11 


2  scalar  and  1  vector  variables  for  each  fluid  and  2  vector  variables  for  the 
electromagnetic  fields. 

The  system  of  equations  can  be  divided  into  the  plasma  fluid  equations  and 
the  electromagnetic  field  equations.  The  plasma  fluid  equations  for  species  a 
are  written  as 
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The  electromagnetic  field  equations  are  written  as 
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Similar  forms  exist  if  either  the  purely  hyperbolic  or  mixed  potential  formula¬ 
tions  are  used. 

The  Jacobians  of  the  hyperbolic  fluxes  dF/dQ  of  the  governing  equations 
[Eqs.  (21)  and  (22)]  are  constructed  in  the  usual  way.  The  eigenvalues  of  the  flux 
Jacobians  give  the  characteristic  velocities.  In  one  dimension,  the  eigenvalues 
of  the  fluid  equations  are 


fluid.  {vaxi  vax  i 


(23) 


where  the  acoustic  speed  for  species  a  is  defined  as 


(24) 


7  is  the  ratio  of  specific  heats,  and  the  temperature  is  defined  from  the  scalar 
pressure,  pa  =  naTa. 

The  electron  acoustic  speed  is  larger  than  the  ion  acoustic  speed  for  the  same 
fluid  temperatures  due  to  the  large  ion  to  electron  mass  ratio.  The  electron 
acoustic  speed  can  be  larger  than  the  Alfven  speed  which  is  a  component  of  the 
eigenvalues  of  MHD.  The  Alfven  speed  for  an  ion/electron  plasma  is  defined  as 


va 


B 

yVo  {rrum  +  mene) 


(25) 


where  B  is  the  magnetic  field  and  /io  is  the  permeability  of  free  space  (47rx  10  7). 
The  eigenvalues  of  the  field  equations  are 


Afield.  =  {±c}  •  (26) 

Therefore,  the  eigenvalues  of  the  multi-fluid  plasma  model  are  generally  not 
bounded  by  the  eigenvalues  of  the  MHD  model.  In  general,  the  fastest  times 
that  must  be  resolved  in  the  multi-fluid  plasma  model  are  the  light  transit  time 
A x/c  and  the  electron  plasma  oscillation  time  . 
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2.2.3  Collisional  Effects 


The  fluids  of  the  multi-fluid  plasma  model  interact  primarily  through  the  elec¬ 
tromagnetic  fields  which  produce  long-range  forces.  However,  the  short-range 
collisional  effects  can  have  a  significant  impact  on  the  complete  plasma  behavior 
and  evolution.  Specifically,  the  collisional  effects  can  effectively  thermalize  the 
fluids  such  that  the  complete  plasma  approaches  a  thermodynamic  equilibrium. 

The  collisional  terms  allow  for  transfer  of  momentum  and  energy  between 
the  different  fluids.  [25]  The  terms  appear  as  frictional  effects  in  Eqs.  (3)  and 
(4) ,  Yip  R-a/3  •  The  rate  of  momentum  transfer  from  species  a  to  species  /3  is 
given  by 

RQ/3  =  mana ua/3  (va  -  v 0) .  (27) 

where  va0  is  the  collision  frequency  between  species  a  and  (3.  Momentum 
conservation  requires  Ra(3  =  —  R/3a  and  Raa  =  0. 


2.2.4  Atomic  Reactions 

The  purpose  of  including  the  effect  of  atomic  reactions  into  the  multi-fluid 
plasma  model  is  to  capture  the  time-dependent  ionization,  recombination,  and 
charge  exchange  reactions  that  are  important  in  laboratory  and  transient  plas¬ 
mas.  Atomic  reactions  lead  to  transitions  between  the  fluids  in  the  multi-fluid 
plasma  model.  For  example,  neutral  ionization  depletes  the  neutral  fluid  and 
increases  the  ion  and  electron  fluids.  The  reactions  also  transfer  momenta  and 
energies  between  the  fluids. 

Contributions  from  atomic  reactions  are  identified  by  terms  on  the  right- 
hand  side  of  Eqs.  (2),  (3),  and  (4).  The  terms  are  expressed  below  for  the 
ionization  (ion),  recombination  (rec),  and  charge  exchange  (cx)  reactions  for  a 
hydrogen  plasma  composed  of  neutral  hydrogen  (a  =  n),  ionized  hydrogen  (a  = 
i),  and  electrons  (a  =  e).  Additional  reactions  and  other  plasma  constituents 
are  also  possible.  The  contributions  to  the  species  densities  are 
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The  contributions  to  the  species  momenta  are 
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The  contributions  to  the  species  energies  are 
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The  reaction  rates  are  given  by 
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where  Twn  is  the  plasma  source  from  electron  impact  ionization  of  neutrals, 
Trec  is  the  plasma  sink  due  to  recombination,  and  Tcx  is  the  plasma  exchange 
due  to  charge  exchange.  The  cross  sections  can  depend  on  temperature. 


2.3  High-Order  Discontinuous  Galerkin  Method 

Electromagnetic  forces  are  exerted  on  the  plasma  fluids  through  the  source 
terms  of  Eq.  (20)  and  the  fluid  motion  affects  the  fields  through  the  source 
terms  of  Eq.  (20).  Even  with  accurate  hyperbolic  flux  calculations,  inaccurate 
calculation  of  the  source  terms  can  lead  to  incorrect  results.  Particularly  in 
equilibrium  situations  where  forces  from  electromagnetic  fields  balance  fluid 
pressure  or  convective  forces,  the  contributions  from  the  source  terms  must  be 
accurately  calculated  to  balance  the  divergence  of  the  hyperbolic  fluxes. 

We  have  developed  a  discontinuous  Galerkin  method  [26-28]  to  solve  the 
governing  equations  of  the  multi-fluid  plasma  model  on  a  computational  grid. 
[2, 6-8]  The  discontinuous  Galerkin  method  is  a  finite  element  approach  that 
allows  for  arbitrarily  high  order  basis  functions  to  model  the  variation  of  the 
system  variables.  Source  terms  are  automatically  coupled. 

The  conserved  variables  of  the  multi-fluid  plasma  model  are  modeled  with 
a  set  of  basis  functions,  vk,  which  can  be  any  desired  order.  The  governing 
equations,  expressed  as  Eq.  (20),  are  multiplied  by  each  basis  function  and 
integrated  over  the  mesh  element  volume  Q.  An  integral  equation  is  generated 
for  each  basis  function. 

f  vk^-dV  +  l  ufcF-dS-  [  F  -VvkdV=  f  vkSdV  (38) 

Jn  dt  Jgn  Jn  Jq 

where  the  divergence  theorem  has  been  applied  to  the  second  term.  The  solu¬ 
tion  variables  Q,  flux  F,  and  source  vector  S  are  given  by  the  first  term,  the 
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second  term,  and  the  right  hand  side  of  Eq.  (20).  The  volume  and  surface  inte¬ 
grals  are  replaced  with  Gaussian  quadrature.  The  source  terms  are  projected 
onto  the  basis  functions  and  are,  therefore,  the  same  order  accurate  as  the  so¬ 
lution  variables.  This  satisfies  the  high-  order  accuracy  requirement  to  preserve 
the  equilibrium  balance  between  the  divergence  of  the  flux  and  the  source. 

The  surface  integral  in  Eq.  (38)  uses  hyperbolic  fluxes  computed  with  a  Roe- 
type  approximate  Riemann  solver.  [1,29]  In  this  method  the  overall  solution  is 
built  upon  the  solutions  to  the  Riemann  problem  defined  by  the  discontinuous 
jump  in  the  solution  at  each  element  interface.  The  numerical  flux  for  a  first- 
order  accurate  (in  space)  Roe-type  solver  is  written  in  symmetric  form  as 

Fi+ 1/2  =  2  C^i+i  +  F%)  ~  2  ^ k  (*2*+i  —  *2*)  l^fcl  Tk  (39) 

k 


where  is  the  kth  right  eigenvector,  is  the  kth  eigenvalue,  and  Ik  is  the  kth 
left  eigenvector,  evaluated  at  the  element  interface  (i  +  1/2).  The  values  at  the 
element  interface  are  obtained  by  a  Roe  average  of  the  neighboring  elements. 
The  flux  calculated  as  above  is  normal  to  the  element  interface  which  is  the 
desired  orientation  for  calculating  the  surface  integral. 

For  the  two-dimensional,  second-order  accurate  algorithm,  a  linear  set  of 
basis  functions  are  used. 


{ut}  =  {vo,vx,vy} 


(  X-Xij  y  -  Vij  \ 
1  ’  Aa™/2  ’  Ay/2  j 


(40) 


where  the  center  of  the  mesh  element  is  located  at  (x{j,  ytJ  )  and  extends  Ax  by 
Ay.  The  conserved  variables  Q  are  defined  as 


Q  —  Qo  Y  Qx^x  T  QyVy  (31) 

within  each  mesh  element.  Update  equations  for  the  coefficients  for  each  con¬ 
served  variable  are  found  directly  from  Eq.  (38)  applied  to  each  mesh  element. 

The  temporal  evolution  is  determined  with  a  Runge-Kutta  method.  A  third 
order  TYD  method  has  been  used  successfully. 

Extensions  to  higher  order  involves  increasing  the  order  of  the  basis  func¬ 
tions  and  following  the  evolution  of  the  additional  coefficients.  For  example  the 
third-order  accurate  basis  functions  would  have  the  form 


—  {vq,  VXi  Vy,  VXy,  VXxt  ^yy} •  (42) 

and  each  conserved  variable  Q  is  defined  as 

Q  —  Q 0  Y  Qx^x  Y  QyVy  T  Qxy^xy  Y  Qxx^xx  Y  QyyVyy  (33) 

within  each  mesh  element. 

The  discontinuous  Galerkin  algorithm  has  been  applied  to  the  electromag¬ 
netic  plasma  shock  demonstrating  the  transition  from  gas  dynamic  shocks  to 
the  MHD  shock  [30, 31]  as  the  Larrnor  radius  is  reduced.  Analysis  of  the  data 
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shows  the  differences  caused  by  the  additional  plasma  waves  that  are  captured 
in  the  two-fluid  model  and,  consequently,  in  the  algorithm  developed  here.  [1] 
It  also  illustrates  the  dispersive  nature  of  the  waves  which  makes  capturing  the 
effect  difficult  in  MHD  algorithms.  The  electromagnetic  plasma  shock  serves 
to  validate  the  algorithm  to  published  data  (MHD  limit)  and  analytical  results 
(gas  dynamic  limit).  The  algorithm  has  also  been  applied  to  study  collisionless 
reconnection  and  the  results  are  compared  to  published  results  of  the  GEM 
challenge  problem.  [24]  The  problem  is  difficult  to  model  and  provides  a  rigor¬ 
ous  test  for  the  algorithm  and  benchmarks  to  other  algorithms.  The  evolution 
of  the  reconnected  magnetic  flux  compares  remarkably  well  with  the  published 
data.  [2]  Additional  applications  are  discussed  in  more  detail  below. 


2.4  Time-Integration  Methods 


The  Runge-Kutta  method  described  previously  has  proven  to  be  robust  and 
accurate  for  explicit  calculations.  The  multi- fluid  plasma  model  has  disparate 
characteristic  speeds  and  frequencies.  The  speeds  range  from  the  high  speed  of 
light  and  electron  plasma  frequency  (nanosecond  times)  to  slow  speed  of  bulk 
fluid  motion  or  ion  sound  speed  (second  or  longer  times).  The  short  time  scales 
dictate  a  small  time-step  for  any  explicit  time-integration  method.  Using  the 
high-oixler  discontinuous  Galerkin  method  further  exacerbates  the  small  time- 
step. 

The  fast  time  scales  contained  in  the  multi-fluid  plasma  model  are  identified 
by  the  eigenvalues  discussed  in  Sec.  2.2.  When  phenomena  of  interest  occurs  on 
these  time  scales  then  time  steps  of  this  size  must  be  used.  When  phenomena 
of  interest  occurs  on  longer  time  scales  then  it  is  desirable  to  use  appropri¬ 
ately  large  time  steps.  As  discussed  earlier,  asymptotic  approximations  can  be 
applied  to  the  multi-fluid  plasma  model  to  remove  the  most  severe  time  scale 
constraints;  however,  these  approximations  fundamentally  change  the  physical 
model  and  yield  an  equation  system  that  can  be  more  difficult  to  solve  mathe¬ 
matically.  We  propose  to  instead  investigate  mathematical  methods  that  allow 
time  steps  larger  than  the  fastest  time  scales. 

Fully  implicit  methods  have  been  investigated.  The  governing  equations  are 
written  as 


dQ 

dt 


f(Q)- 


(44) 


The  time  advance  can  be  expressed  in  an  implicit  form  of  arbitrary  order  ac¬ 
curacy  as 


Qn+1  =  g  (At,  Qn+\  Qn,  Qn~\  ...)  .  (45) 


A  second-order  Crank-Nicolson  method  has  been  implemented,  where  the  time 
advance  equation  is  written  as 


Qn+1  =  Qn  +  T[f  {Qn+1)  +  /  (Qn)} , 


(46) 


which  is  solved  iteratively  using  a  variety  of  implicit  numerical  methods,  e.g. 
Conjugate  Residual,  Biconjugate  Gradient,  and  GMRES,  with  preconditioners 
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Figure  3:  Electron  and  ion  number  densities  for  the  electromagnetic  shock  problem 
solved  with  the  implicit  method  using  time  steps  200  times  larger  than  the  fastest 
time  system  time  scale.  An  explicit  solution  is  also  shown  for  comparison. 

such  as  Incomplete  LU  and  Additive  Schwarz.  The  method  has  been  applied 
to  several  benchmark  problems  to  compare  the  accuracy  of  the  solution  when 
large  time  steps  are  used.  Solution  accuracy  is  well  preserved  for  the  resolved 
timescale.  However,  the  full  implicit  time  integration  is  inefficient  because  of 
the  large  size  of  the  equation  system. 

The  electromagnetic  shock  problem  provides  a  rigorous  benchmark  because 
it  has  fast  phenomena  associated  with  light  and  whistler  wave  propagation, 
and  it  has  slower  phenomena  associated  with  the  massive  ion  motion.  Figure  3 
shows  simulation  results  for  the  ion  and  electron  number  densities  for  the  elec¬ 
tromagnetic  shock  problem  where  steps  200  times  larger  than  the  light  transit 
time  and  20  times  larger  than  the  electron  characteristic  time  are  used.  An 
explicit  solution  which  resolves  all  of  the  time  scales  is  also  shown  for  compar¬ 
ison.  The  solution  shows  that  the  electron  dynamics  are  only  approximately 
modeled  by  the  implicit  method  as  expected,  but  the  ion  dynamics  are  well  cap¬ 
tured  and  agree  with  the  explicit  solution.  Both  methods  use  a  second-order 
discontinuous  Galerkin  spatial  representation. 

A  semi-implicit  method  has  been  formulated  to  perform  a  splitting  of  the 
governing  equations  of  the  two-fluid  plasma  model.  The  splitting  is  based  on 
the  different  expected  physics  and  corresponding  different  scales  dictated  by 
the  specific  model  equations.  Specifically,  the  electron  fluid  and  electromag¬ 
netic  fields  introduce  the  fast/small  scales  and  are  solved  implicitly  using  the 
method  described  by  Eq.  (46),  and  the  ions  introduce  slow/large  scales  and 
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are  solved  explicitly  using  the  third-order  accurate  TVD  Runge-Kutta  time  ad¬ 
vance  method.  The  combination  of  implicit  and  explicit  treatments  limits  the 
size  of  the  operator  matrix  while  still  allowing  time  steps  that  are  not  limited 
by  the  fast  timescales.  The  time  step  is  still  limited  by  the  ion  motion;  however, 
the  ion  motion  is  often  the  timescale  of  interest  and  needs  to  be  resolved  for 
accurate  simulation  results.  The  semi-implicit  method  has  yielded  encouraging 
results  by  exploiting  the  fundamental  physics  of  the  governing  equations.  How¬ 
ever,  there  are  additional  avenues  to  pursue  to  further  incorporate  the  physics 
of  the  governing  equations  into  the  numerical  method. 

2.5  Applications 

The  described  above  is  implemented  on  parallel  computers  using  an  automatic 
domain  decomposition  technique  with  MPI  message  passing.  The  algorithm 
is  implemented  in  a  code  called  WARPX  (Washington  Approximate  Riemann 
Plasma)  which  uses  C++  object  oriented  programming  and  other  modern  soft¬ 
ware  techniques  to  simplify  the  maintainability  and  extendibility  of  the  code 
and  HDF5  for  parallel  output. 

2.5.1  Field  Reversed  Configurations  (FRC)  in  Three  Dimen¬ 
sions 

WARPX  has  been  applied  to  study  Z-pinch  dynamics  and  three-dimensional 
FRC  evolution  to  investigate  anomalous  resistivity  that  experimentally  limits 
the  plasma  current.  An  experiment  at  Kirtland  AFB  (AFRL/RDHP)  forms 
FRC  plasmas,  translates  them  into  a  cylindrical  flux  conserver,  and  compresses 
the  plasma  to  high  energy  densities  by  imploding  the  flux  conserver  onto  the 
plasma.  FRC  plasmas  have  been  modeled  using  MHD  codes  with  unsatisfac¬ 
tory  results.  Specifically,  FRC  plasmas  are  characterized  by  large  gyroradii 
(. L/tl  =  2  —  4)  and  charge  separation.  MHD  models  have  failed  to  predict 
many  experimental  observations,  such  as  observed  stability  and  anomalous  re¬ 
sistivity.  The  lack  of  agreement  is  conjectured  to  be  caused  by  two-fluid  effects. 
Sample  evolutions  of  the  FRC  plasmoid  are  shown  in  Fig.  4. 

2.5.2  Three-Fluid  Plasma  Model  for  Plasma  Production 

WARPX  provides  a  flexible  code  framework  that  allows  easy  extension  of  the 
physical  model  to  include  multiple  fluids.  We  have  extended  WARPX  to  a 
three-fluid  (electrons,  ions,  and  neutrals)  simulation  of  plasma  sheath  forma¬ 
tion.  Atomic  reactions  are  incorporated  that  describe  the  effects  of  collisions 
between  the  species  explicitly,  allowing  for  the  identification  of  regions  of  ioniza¬ 
tion/recombination  and  interspecies  momentum  and  energy  transfer.  Plasma 
sheath  formation  is  important  for  electrode-based  plasma  technologies,  e.g. 
plasma  actuators  for  control  of  high-speed  aerospace  vehicles.  The  multi-fluid 
model  captures  electron  inertial  effects  and  has  revealed  a  new  physical  effect. 
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Figure  4:  Three-dimensional  evolution  of  an  FRC  using  the  WARPX  code.  Small 
scale  variations  are  evident  in  the  azimuthal  direction  indicative  of  a  drift  instability. 

During  the  initial  formation  of  the  plasma  sheath  the  applied  electrode  po¬ 
tential  excites  a  Langmuir  wave  that  propagates  into  the  bulk  plasma.  The 
dispersion  relation  is  given  by 

w2  =  Wpe  +  \k2v^e.  (47) 

Typically,  sheath  simulations  assume  electrostatic  fields  and  miss  the  elec¬ 
trodynamics  of  the  formation  process.  Propagating  Langmuir  waves  are  shown 
in  simulation  results  in  Fig.  5.  The  numerical  dispersion  agrees  with  the  dis¬ 
persion  relation  of  Eq.  (47). 

The  plasma  sheath  formation  studies  provide  better  physical  understanding 
into  the  plasma  production  process.  Plasma  (ions  and  electrons)  is  produced 
by  ionizing  the  neutral  gas,  and  it  is  lost  when  it  reaches  the  electrode  and 
recombines.  An  analytical  model  of  the  electrode  describes  secondary  electron 
emission  and  recombination  at  the  electrode.  In  addition  to  the  plasma  sheath 
that  naturally  forms  around  electrodes,  a  voltage  can  be  applied  to  the  elec¬ 
trodes  to  drive  a  current  through  the  plasma.  Phenomena  such  as  cathode  and 
anode  drops  are  accurately  simulated. 

3  Conclusions 

Investigating  advanced  plasma  models  is  motivated  by  the  need  to  simulate 
complicated  plasma  physics  phenomena  that  is  not  captured  in  simpler  models. 
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Figure  5:  Initial  plasma  production  showing  the  propagation  of  a  Langmuir  wave  in 
the  electron  density  (black  trace)  and  electric  field.  Numerical  dispersion  calculated 
in  (k,oj)  space  agrees  with  analytical  dispersion  relation  (dashed  line). 
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The  multi-fluid  plasma  model  is  proving  to  be  a  model  that  is  significantly  more 
advanced  and  complete  than  the  usual  MHD  model  without  the  computational 
complexity  required  in  general  kinetic  models.  This  project  has  developed  the 
multi-fluid  plasma  model  and  cast  the  governing  equations  in  a  balance  law  form 
that  lends  itself  to  accurate  numerical  solutions.  The  algorithm  developed  in 
the  project,  its  implementation  into  WARPX,  and  its  application  to  benchmark 
and  real  experimental  problems  have  demonstrated  the  capability  of  both  the 
multi-fluid  plasma  model  and  the  numerical  techniques  in  the  algorithm. 
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